clc;
clear;

%% main主脚本
tic;
startTime = cputime;
pde = process_pde_elliptic(1);
process_plot_ex_scalar(pde, pde.u);
vector_h = [1/2 1/4 1/8 1/16];
basis_type = "P1";

u_Linf = zeros(size(vector_h,2),1);
u_L2 = zeros(size(vector_h,2),1);
u_H1 = zeros(size(vector_h,2),1);

for i = 1:size(vector_h,2)
    h = vector_h(i);
    mesh = process_mesh(pde, h);
    fs = process_fespace(mesh, basis_type);

    % process_showmesh(mesh, fs, [1 1 0 0]);

    A = term_grad_P_grad_Q(mesh, fs);
    b = term_F_Psi(mesh, fs, pde.f);
    [A, b] = process_boundary_elliptic(A, b, mesh, fs, pde);
    SOL = A\b;

    u_Linf(i) = process_error_elliptic(mesh, fs, pde, SOL, "Linf");
    u_L2(i) = process_error_elliptic(mesh, fs, pde, SOL, "L2");
    u_H1(i) = process_error_elliptic(mesh, fs, pde, SOL, "H1");

    fprintf("已完成\t%d/%d\n", i, size(vector_h,2));
end
process_plot_scalar(mesh, fs, SOL);

process_rate(u_Linf, vector_h);
process_rate(u_L2, vector_h);
process_rate(u_H1, vector_h);

TOCtime = toc;
CPUtime = cputime - startTime;
disp(["TOC time(s)", num2str(TOCtime)]);
disp(["CPU time(s)", num2str(CPUtime)]);


%% process_pde_elliptic
function pde = process_pde_elliptic(index)
switch index
    case 1
        pde.index = index;
        pde.u = @(x,y,z) cos(pi*x).*cos(pi*y).*cos(pi*z);
        pde.udx = @(x,y,z) -pi*sin(pi*x).*cos(pi*y).*cos(pi*z);
        pde.udy = @(x,y,z) -pi*cos(pi*x).*sin(pi*y).*cos(pi*z);
        pde.udz = @(x,y,z) -pi*cos(pi*x).*cos(pi*y).*sin(pi*z);
        pde.f = @(x,y,z) 3*pi^2*cos(pi*x).*cos(pi*y).*cos(pi*z);
        pde.domain = [0, 1, 0, 1, 0, 1];
        pde.bdary = [0, 1, 0, 0, 0, 0];
    case 2
        pde.index = index;
        pde.u = @(x,y,z) exp(x+y+z);
        pde.udx = @(x,y,z) exp(x+y+z);
        pde.udy = @(x,y,z) exp(x+y+z);
        pde.udz = @(x,y,z) exp(x+y+z);
        pde.f = @(x,y,z) -3*exp(x+y+z);
        pde.domain = [-0.5, 0.5, -0.5, 0.5, -0.5, 0.5];
        pde.bdary = [0, 0, 0, 0, 0, 0];
    case 3
        pde.index = index;
        pde.u = @(x,y,z) x.*y.*z.*(1-0.5*x).*(1-y).*(1+0.5*z).*exp(x+y+z);
        pde.udx = @(x,y,z) (y.*z.*exp(x+y+z).*(x.^2-2.0).*(y-1.0).*(z+2.0))./4.0;
        pde.udy = @(x,y,z) (x.*z.*exp(x+y+z).*(x-2.0).*(z+2.0).*(y+y.^2-1.0))./4.0;
        pde.udz = @(x,y,z) (x.*y.*exp(x+y+z).*(x-2.0).*(y-1.0).*(z.*4.0+z.^2+2.0))./4.0;
        pde.f = @(x,y,z) (y.*exp(x+y+z).*(x.*-1.2e+1-z.*4.0-x.^2.*z.^2+x.*y.*1.2e+1+x.*z.*4.0+y.*z.*4.0-x.^2.*y.*6.0+x.*z.^2.*6.0+x.^2.*z.*2.0+y.*z.^2.*2.0+x.^2.*6.0-z.^2.*2.0+x.*y.*z.^2.*2.0-x.^2.*y.*z.*1.0e+1-x.^2.*y.*z.^2.*3.0+x.*y.*z.*1.2e+1))./4.0;
        pde.domain = [1, 2, 0, 1, -2, -1];
        pde.bdary = [0, 0, 0, 0, 0, 0];
    otherwise
        error("Invalid pde index.");
end
fprintf("pde index: %d\n", pde.index);
fprintf("pde domain: xmin,xmax,ymin,ymax,zmin,zmax\n");
disp(pde.domain);
fprintf("pde boundary: left,front,right,back,bottom,top (0=Dirichlet,1=Neumann)\n");
disp(pde.bdary);
end

%% process_boundary_elliptic
function [A, b] = process_boundary_elliptic(A, b, mesh, fespace, pde)
Nb = fespace.Nb;
Nlb = fespace.Nlb;
% Neumann
for k = 1:size(pde.bdary,2)
    if 1 == pde.bdary(1,k)
        bdf = bdary_face(mesh, k);
        rows = zeros(size(bdf,1)*Nlb,1);
        val = rows;
        flag = 1;
        for n = 1:size(bdf,1)
            points = mesh.node(bdf(n,2:4),:);
            AB = points(2,:)-points(1,:);
            AC = points(3,:)-points(1,:);
            vector_n = cross(AB, AC);
            vector_n = -vector_n./norm(vector_n);
            fb = @(x,y,z) vector_n(1)*pde.udx(x,y,z) + vector_n(2)*pde.udy(x,y,z) + vector_n(3)*pde.udz(x,y,z);
            e = bdf(n,1);
            V = mesh.node(mesh.elem(e,:),:);
            [int_points_loc, int_weights_loc] = integral_loc(points, "2d");
            phi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace.basis_type, "x");
            ffb = fb(int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3));
            for i = 1:fespace.Nlb
                f = ffb .* phi(:,i);
                r = int_weights_loc * f;
                rows(flag) = fespace.elem(e,i);
                val(flag) = r;
                flag = flag + 1;
            end
        end
        vector_b = sparse(rows, 1, val, Nb, 1);
        b = b + vector_b;
    end
end
% Dirichlet
for k = 1:size(pde.bdary,2)
    if 0 == pde.bdary(1,k)
        bdn = bdary_node(mesh, fespace, k);
        A(bdn,:) = 0;
        A(sub2ind(size(A), bdn, bdn)) = 1;
        b(bdn) = pde.u(fespace.node(bdn,1),fespace.node(bdn,2),fespace.node(bdn,3));
    end
end
end

%% process_error_elliptic
function error_u = process_error_elliptic(mesh, fespace, pde, U, error_type)
switch error_type
    case "Linf"
        err_u = 0;
        for e = 1:mesh.Ne
            V = mesh.node(mesh.elem(e,:),:);
            [int_points_loc, ~] = integral_loc(V, "3d");
            phi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace.basis_type, "x");
            uh = phi * U(fespace.elem(e,:));
            f = max(abs(pde.u(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - uh));
            if f > err_u
                err_u = f;
            end
        end
        error_u = err_u;
    case "L2"
        err_u = 0;
        for e = 1:mesh.Ne
            V = mesh.node(mesh.elem(e,:),:);
            [int_points_loc, int_weights_loc] = integral_loc(V, "3d");
            phi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace.basis_type, "x");
            uh = phi * U(fespace.elem(e,:));
            f = (pde.u(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - uh).^2;
            r = int_weights_loc * f;
            err_u = err_u + r;
        end
        error_u = sqrt(err_u);
    case "H1"
        err_ux = 0;
        err_uy = 0;
        err_uz = 0;
        for e = 1:mesh.Ne
            V = mesh.node(mesh.elem(e,:),:);
            [int_points_loc, int_weights_loc] = integral_loc(V, "3d");
            phix = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace.basis_type, "dx");
            phiy = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace.basis_type, "dy");
            phiz = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace.basis_type, "dz");
            uhx = phix * U(fespace.elem(e,:));
            uhy = phiy * U(fespace.elem(e,:));
            uhz = phiz * U(fespace.elem(e,:));
            fx = (pde.udx(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - uhx).^2;
            fy = (pde.udy(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - uhy).^2;
            fz = (pde.udz(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - uhz).^2;
            rx = int_weights_loc * fx;
            ry = int_weights_loc * fy;
            rz = int_weights_loc * fz;
            err_ux = err_ux + rx;
            err_uy = err_uy + ry;
            err_uz = err_uz + rz;
        end
        error_u = sqrt(err_ux + err_uy + err_uz);
    otherwise
        error("Invalid error type.");
end
end